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Abstract 

A finite-difference machine code is brought to bear on the 
wake-vortex problem in the quasi— cylindrical boundary-layer approxi' 
nation. A turbulent-energy model containing new features is 
developed that accounts for the major effects disclosed by more 
advanced models in which the parameters are not yet established. 
Several puzzles are resolved that arose in previous theoretical 
investigations of wake vortices. 
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Introduction 

Methods for reduction of velocities in trailing vortices of 
large aircraft are of current interest for the purpose of shortening 
the waiting time between landings at central airports. Procedures 
for estimating the response of the vortices to stimuli such as air 
injection or special flap deflections are required. The favorable 
effect of a special ogee extension at the wing tip has been investi- 
gated in wind tunnel experiments by Rorke and Moffitt. It is 
necessary to determine the effects of this and possibly numerous j 
other proposals in the far flow field as well as the near field 
that can be investigated in wind tunnels. Computer simulation of 
this problem may eventually produce competitive flexibility. 

At the present time definitive simulation of the rate of decay 

and growth in core size due to turbulent viscosity has not been 

achieved. A useful review of previous work on this problem and 

contributions to its solution have been given by Govindaraju and 

Saf fman. ^ Several attempts have been made to stimulate the effects 

of turbulence in wake vortices using integral methods, notably 
3 4 

Kuhn and Nielsen and Fernandez and Lubard. Additional useful 
information resulted from these efforts. The wind-tunnel measure- 
ments of Hoffman and Joubert 5 and Chigier and Corsiglia 6 have been 
used extensively and have contributed to the present knowledge of 
turbulent wake vortices. A preliminary calculation of the flow 
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in an isolated vortex has been made by Donaldson^ based on a sophis- 
ticated turbulence theory. The method of invariant modeling employed 
should be useful in the development of future turbulence theories. 

Investigators at the Imperial College in London (e.g., Launder 
et al.^) have developed turbulence models in which the length scale 
is governed by an additional differential equation in lieu of assign- 
ment of a length scale as In Donaldson's model. Application of this 
concept to turbulent wake vortices in the present paper leads to 
significant departure from Donaldson's prediction of the rate of 
decay of such vortices. An effect of suppression of turbulence 
due to curvature of the flow contained in Donaldson's model) but 
not in that | of Launder et al. , also plays an important role. A 
method for including these two effects in an extension of Prandtl's 
turbulent energy model is developed in the present paper. The 
relatively slow rate of decay of turbulent wake vortices that has 
been observed in flight experiments 10 is thereby explained. The 
approach to self-similar solutions predicted in earlier investiga- 
tions by Donaldson and Baldwin et al. does not occur when both 
of the above effects are considered. 

Notation 

a^ parameter in turbulent energy model (0.15) 

C parameter (0.5 C ,C = 0.0643) 

X t j. y 

C 2 parameter = 0.223) 


- 3 - 




max 



2 


q 


! 


r,z 


r 


1 


r 


2 


R,Z 

s 

S 


u,v,w 


parameter (1.0) 

parameter (1.43) 

parameter (1.92) 

2 

parameter (4a^ = 0.09) 

dimensionless kinematic viscosity (Eq. (7)) 
maximum value of f in a given profile 
mixing length 

Prandtl number in dissipation equation (= 1.3) 

Prandtl number in energy equation (= 1.0) 
dimensionless pressure (Eq. (3)) 

approximate dimensionless total pressure (Eq. (6)) 
specific (amount per unit mass) time— averaged turbulent 
fluctuation kinetic energy (times 2) 
dimensionless radial and axial coordinates 
(r = W R/r , z = W z/r ), (e.g. Eq. (l)-(5)) 
dimensionless radius at the position of maximum circumfer- 
ential velocity (e.g. Eqs. (16) and (17)) 
dimensionless radius at outer boundary of computational 
mesh (e.g. Eq. (9)) 

dimensional radial and axial coordinates 
dimensionless rate of strain 
dimensional rate of strain 

dimensionless radial, circumferential, and axial velocity 
components (e.g., u = U/w^), (e.g., Eqs. (1) - (5)) 
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U,v,w 


dimensional radial, circumferential and axial velocity 
components 

W dimensional free stream velocity 

00 

r circulation (= 2-rrVR) 

circulation at position of maximum circumferential velocity 

r circulation at large R 

00 

Ar computational mesh size of r coordinate (- 

Az variable computational step size 

e specific dissipation of turbulent energy 

v kinematic viscosity 

v turbulent kinematic viscosity 

T 

Theory 


Quasi-Cylindrical Flow 

The basic equations in the boundary-layer approximation are 

taken from Lilley and Chigier. 12 Attention has been confined to 

incompressible flow in this paper and the equations reduced to a 

nondimensional form through division of velocities by the free-stream 

2 

velocity W^, division of pressure by , and multiplication of 

axial and radial coordinates by the factor W^/T^. The resulting 
equations in a nondimensional form suitable for finite difference 
solution are 



-u 


9w 1 f 9w 
9r r 3r 3r 


( 1 ) 
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where 



2 

w 



( 6 ) 

(7) 


The turbulent viscosity coefficient is determined from a 

turbulence model to be discussed. The boundary condition are 



u ( r 2 ’ z ) = 0 

(8) 


v ( r 2 ,z ) = 2»r 2 

(9) 


W (r^,z) = 1 

(10) 


P ( r 2’ z ) = 0 

(10a) 

where 

r^ is a radius large compared to the size of the core 

(r x ). 

At the 

center of the core 
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u(0,z) = 0 
v(0,z) = 0 

The initial conditions are 


u ( r ’ z 0 ) = 0 



(ID 

( 12 ) 


(13) 

(14) 


(15) 


where r^ is the radius at which the circumferential velocity is a 
maximum. 

Turbulence Model 

Since detailed experimental information on the flow in turbulent 
vortices is not yet available, methods that depend on evaluation of 
parameters from the particular experiment at hand cannot be used. 
Instead a complete model is needed that bears an approximate relation- 
ship to the exact flow in a manner analogous, for example, to the 
relationship between simple kinetic theory and the rigorous kinetic 

g 

theory of gases* It has been shown by Launder, et al. that Prandtl's 

9 

1945 turbulent energy model works reasonably well in axisymmetric 
wakes. The basic equation can be written (in the notation of A. 

Towns end 
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Dt 


PR, 


$ * V 


V) + y^ 2 ' 


ST 3/2 /-j\ 3/2 

q - a, (q J /L + F 


(16) 


where F is an additional term developed in this paper. The 
kinematic-turbulent-viscosity coefficient and local shear are 

given by 



(17) 

(18) 


For later discussion it is worthwhile to consider the production 
term of equation (16) in more detail. Since 

2 2 y 2 

R 2 R' 


LL V \ 2 _ / 9 V\ 2 + 3 v! _ 2 X 
\3R R j \3R/ r 2 2 


m 

3R 


(19) 


The production term can be written 

2 


V^ Lsi 


r 2 

q = 


/3W\ 2 /3V\ 2 , _ V 2 V, 3 (VR) r~2 

\Ti) + \W +3 ^2 - 2 ^2^r~ L V q 


In axisymmetric wakes, the mixing length L is assigned a uniform 
value that is taken to be a constant fraction of the width of the 
wake. However, Launder, et al., have developed an additional 
equation from which the mixing length can be computed rather than 
assigning its value (again in Townsend's notation) 


■ PT * * V T^ + c i^ 2s2 - C 2V < 2e/L 


( 20 ) 


( 21 ) 
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and 


L=a 3/2 ( q 2) 3/2 /e 


( 22 ) 


8 14 

Launder, et al. and Hanjalic and Launder have shown that computation 
of the mixing length in this manner leads to Improvement of comparisons 
with measurements in a variety of turbulent free shear flows as well 
as boundary layers. 

The foregoing turbulence model does not apply to turbulent 
vortices because it does not account for suppression of turbulence 
due to curvature of the mean flow that has been disclosed by more 
advanced models such as that of Donaldson^. This effect is illustrated 
in figure 1, which represents a rotating flow in cross section. The 
Rayleigh stability criterion indicates that if the product VR decreases 
in the outward direction, the flow is unstable against disturbances 
that would interchange fluid between inner and outer regions. Con- 
versely, if the product VR increases in the outward direction, there 
is a damping effect on turbulent eddies produced by the shearing motion. 
This concept can be made quantitative by finding the energy per unit 
mass of fluid required to interchange fluid between the inner and 
outer annuli, while retaining an unchanged angular momentum of the 
fluid transferred 


12 1 2 
A^q = 2 ^ (VR) * 



, V 3VR 2 
r2 3L 
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The rate of transfer of fluid by an eddy of size L and velocity 

~ 

q is approximated by 


1_ dm „ . 
m dt ~ \ 



These relations can be combined to obtain an additional term to be 
included in the turbulent energy equation (18) ♦ The new turbulence 
suppression term is 

F ■ - 2 % Hr 1 L V^ c 3 <23> 

This relation applies more generally to flows with curvature if R 
is taken to be the local radius of curvature of streamlines. Present 
information does not permit determination of a better value of the 
magnitude of this term than that resulting from the above derivation 

(c 3 = 1.0). 

It is interesting to note that the suppression term is of the 
same form and same sign as a term arising from the shear that is 
already present in Prandtl r s energy equation, namely the last term 
in eq. (20) • One might wonder whether the curvature suppression 
effect discussed above is already included in Prandtl r s model. How- 
ever, as can be seen in equation (18), the shear contribution to the 
production of turbulence must always be positive because it consists 
of a sum of squared quantities. It seems clear that Rayleigh in- 
stability should enhance the growth of turbulence and an opposite 
stabilizing condition of the mean flow should be capable of stopping 



the growth entirely as in the suppression of turbulence due to 
buoyancy forces (e.g., Townsend^) - We therefore conclude that the 
term given in equation ( 23 ) should be included again even if it does 
already appear in the production term due to a different mechanism 
(shearing motion rather than curvature of the mean flow) . 

It may be possible to simulate the suppression of turbulence 
by making the length L smaller so that the effect of turbulent 
mixing is reduced and the rate of production due to shear decreased. 
An improvement in the agreement of calculations with measurements in 
curved boundary layers was achieved in this manner by Bradshaw^. 

In near-equilibrium situations it may be possible to relate that 
method to the present one. However, the mechanism by which such an 
effect would enter the dissipation equation (eq. ( 19 )) is not clear 
to the present authors, whereas the origin of the energy suppression 
term (eq. ( 23 )) seems relatively transparent. 

As mentioned earlier, the turbulence model of Donaldson^ con- 
tains a curvature suppression effect that enters through the tensor 
formalism used in the method of invariant modeling. Ultimately, such 
an approach may be more satisfactory than that employed here. We 
have been unable to relate our method to that of Donaldson and have 
not been convinced that his treatment of the effect is yet in final 
form. For the present, we prefer simple physical reasoning applied 
to the special type of flow under consideration in which the origin 
of the effect seems understandable. 
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The nond linens ional forms of the energy and dissipation equations 


used in the finite-difference calculations are 


where 


_3K _ _1 
3z w 


3K 1 X 

U 3r PR q r 3r 
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The dimensionless local shear s and length l are given by 


s = 




l 


3/2 

K 

K 


(28) 

(29) 


g 

The numerical values of the constants (taken from Launder et al. ) 
are given in the list of notation. The dimensionless turbulent 
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viscosity coefficient is 



r 

OP 


= C Ski 

y 


( 30 ) 


Computational Procedure 

A machine code was developed for finite-difference solution of 

the foregoing relations according to the explicit predictor-corrector 

method of MacCormack . ^ Equations (1), (2), (24), and (25) are stepped 

forward in z to obtain values of p , v, k., and £. At each step 

w 

Eq. (3) is integrated by quadrature to obtain p. Equation (4) can 
then be used to evaluate w. The derivative 3w/9z\ is computed by 
finite-difference and Eq. (5) Integrated by quadrature to obtain u. 

In the explicit method used to step forward in z all auxiliary 
quantities are evaluated In terms of the Independent and primary de- 
pendent variables by means of the remaining formulas. One hundred 
equally spaced mesh points are used for the radial coordinate r. The 
step size of the axial coordinate z is dictated by the requirements 
of numerical stability. Instabilities arising from the last term 

in Eq. (1) are dominant according to an analysis based on the methods 

18 

of Lomax and Bailey. A workable stability criterion was found to 
be 

2 

Az < 0.3(Ar) ^2) 

^max 
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RESULTS 


For isotropic shear flows the simplified turbulence model described 
accounts for the major effects disclosed by more advanced models in 
which the parameters are not yet established. It therefore applies 
to turbulent vortices with small axial velocities and can be expected 
to lead to realistic predictions of changes taking place in the far 
flow field of wake vortices. We have made finite difference calcu- 
lations based on this model in the quasicylindrical approximation. 
Information on initial values of the axial and radial velocities 
being unavailable, they have been set equal to zero. Thus calculations 
in this report correspond to an isolated vortex extending from to 
+«> through the transformation t = Z/W^. 

In the middle of figure 2 are shown computed circulation profiles 
with free-stream velocity W^, circulation at large radius F^, and 
wing chord c assigned values corresponding to the Cherokee flight 
data of McCormick, Tangier, and Sherrieb^. The calculations were 
started at z = 0 with an initial circulation profile corresponding 
to that of a self-similar laminar vortex (Lamb vortex) . As shown in 
the profile at 34.5 chord lengths behind the wing, an overshoot 
develops in the circulation and moves outward. In the inner regions 
where the circulation rises steeply the curvature suppression results 
in reduced mixing so that changes take place more slowly than in the 
outer regions. The corresponding velocity profiles appear at the top 
in figure 2. 
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At the bottom of figure 2 a typical turbulent intensity profile 
at a downstream station 34.5 chord lengths behind the wing is shown. 
The solid curve represents the computed turbulent intensity and the 
dashed curve indicates the equilibrium values that would be present 
if the production and dissipation terms in the turbulent energy 
equation (eq. (18))^ were balanced to zero. The equilibrium values 
are large at a radius of 40 inches because the gradient of circula- 
tion is negative in that region as can be seen in the middle curves 
of figure 2. This negative gradient leads to an enhancement rather 
than suppression of turbulence as indicated by equations (18) J and (23). 
However, there is a lag of the actual turbulence below the equilibrium 
values. The turbulence is weak inside the radius of maximum velocity 
so that relatively slow changes in velocity occur there. The local 
peak in turbulent intensity near the radius of maximum velocity can 
be attributed to changes outside this radius which enhance the local 
shear leading to production of turbulence. 

Figure 3 contains profiles of the computed mixing length at 
three axial stations. The calculation was started with a uniform 
mixing length at z - 0. However, the values quickly change within 
10 chord lengths and become highly nonuniform across the profile. 

In the region of circulation overshoot where the turbulence is en- 
hanced, the eddies become quite large. The eddy size quickly goes 
to very small values inside the radius of maximum velocity where the 
turbulent intensity is also small due to curvature suppression. In 
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figure 4 are shown turbulent viscosity profiles. These results 
indicate that there is considerable turbulent mixing in the outer 
part of the vortex, but very little near the point of maximum 
velocity. 

Figure 5 contains estimates of the decay of maximum circumferen- 
tial velocity according to several turbulence models. The circles 

represent the flight data of McCormick et al.^. The lower curve 

4 11 

indicates a typical result from previous investigations * in which 
the suppression of turbulence due to curvature of the flow is not 
accounted for. In the middle curve, this effect was included but 
the mixing length was taken to be a constant fraction of the radius 
of maximum velocity. In that case the calculation indicates appre- 
ciable encroachment of turbulence into the core by diffusion as in 
Donaldson's calculation^. However, if the variable mixing length 
is accounted for by integration of equation (21) , such diffusion is 
inhibited by the very small mixing length that exists near the radius 
of maximum velocity. The top curve is based on the complete model 
with all large effects accounted for. 

There are several possible explanations for the discrepancy 
between the upper curve and the flight data represented by the 
circles. Calculations based on different initial circulation pro- 
files may lead to better agreement. Also it is known that the level 
of ambient atmospheric turbulence affects the rate of decay. Rapid 
decays occur in rough air possibly due to time dependent distortions 
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of the circulation profiles introduced by the buffeting. In calm 
air slow decays take place that would be in better agreement with 
calculations not including such effects. The important point is 
that the present simplified turbulence model predicts slow decays 
in contrast to previous results that were in gross disagreement with 
flight data. This model indicates that the level of turbulence 
inside the radius of maximum velocity is much lower than previously 
supposed. The initial roll up is more nearly inviscid than many of 
us believed, particularly in the inner regions. Two basic mechanisms 
are responsible for this. One is the suppression of turbulence due 
to curvature of the flow (eq. (23)). The other is the variable 
mixing length (fig. 3) which becomes small near the radius of maximum 
velocity, thereby inhibiting the inward diffusion of turbulence. 

These effects lead to small turbulent mixing in the inner regions of 
the vortex such that the approach to self-similar solutions noted 
in previous investigations^*^ does not occur. 
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Figure Titles 


Fig, 1 Rotating flow in cross section. 

2 

Fig. 2 Profiles in wake vortex (W = 132 ft/sec, r = 112 ft /sec, 

v OO CO 

c = 5.26 ft). 

Fig. 3 Mixing length (W ro = 132 ft/sec, = 112 ft^/sec, c = 5.26 ft). 

2 

Fig. 4 Turbulent viscosity profiles (W^ = 132 ft/sec, = 112 ft /sec, 
c - 5.26 ft). 

2 

Fig. 5 Decay of wake vortex (W^ = 132 ft /sec, = 112 ft / sec, 
c = 5.26 ft) . 
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